3 Annotating land use

In this study we aim to quantify the response to fireworks across different species groups, for which take-off habitat is probably a good indicator. In this notebook we will classify land use, and a variety of factors that can influence the ‘intensity’ of fireworks disturbance, such as the distance to the nearest urbanised area for each of the PPI ‘pixels’.

The land use is based on the CORINE Land Cover dataset and specifically the 2018 version (CLC2018), which should be most relevant for the 2017-2018 fireworks event.

3.2 Converting the land use map

To start, we need to convert the land use map to the same 1) resolution, and 2) extent of the radar PPIs as we can then simply ‘overlay’ both rasters on top of each other and do calculations. We load the PPIs and extract the CRS information contained in the proj4 strings.

And we load and prepare the land use map it’s all about. To aid the classification process, we also load the land use classes contained in the entire CLC2018 dataset, otherwise the classes will remain anonymous numbers.

3.2.1 Cropping the land use map

As the CLC2018 dataset is so large it does not fit in memory at all in the steps below, so we have to crop the raster dataset for further processing. Even then, it still requires a beefy computer to process these files. First we calculate a bounding box for the landuse raster based on the bounding boxes of the radar data.

We can now crop and plot the land use maps centered on the radar sites, with a 2.510^{4} meter padding surrounding the extent of the radar data.

And plot the result:

3.2.2 Reprojecting the land use map

Now that the land use map is cropped, we can start the reprojection to the CRS of the radar PPI. As we’re dealing with categorical data, we set method = "ngb" to use nearest neighbour interpolation. Despite using interpolation, this reprojection does not change the resolution from the base of the CLC2018 dataset to that of the PPI, so we’ll have to do that next.

If the reprojection went successful, the CRS of the reprojected land use map and the radar PPI should be the same.

## [1] TRUE
## [1] TRUE

Apparently that is the case.

3.2.3 Reclassifying land use types to functional classes

The CLC2018 dataset contains a total of 44 land use classes. For our purpose, we reduce these 44 to 5 classes with more biologically relevant groupings, specifically:

  1. Urban area
  2. Agricultural area
  3. Semi-open area
  4. Forests
  5. Wetlands
  6. Water bodies

The following table is used to convert/reclassify the classes contained within the CLC2018 dataset, with the original land use classes under landuse_class and how they will be reclassified under landuse_target. We also indicate whether these areas are inhabited (inhabited = 1), or uninhabited (inhabited = 0).

landuse_id landuse_class landuse_target landuse_target_id inhabited
111 Continuous urban fabric Urban area 1 0
112 Discontinuous urban fabric Urban area 1 1
121 Industrial or commercial units Urban area 1 0
122 Road and rail networks and associated land Urban area 1 0
123 Port areas Urban area 1 0
124 Airports Urban area 1 0
131 Mineral extraction sites Urban area 1 0
132 Dump sites Urban area 1 0
133 Construction sites Urban area 1 0
141 Green urban areas Urban area 1 0
142 Sport and leisure facilities Urban area 1 0
211 Non-irrigated arable land Agricultural area 2 0
212 Permanently irrigated land Agricultural area 2 0
213 Rice fields Agricultural area 2 0
221 Vineyards Semi-open area 3 0
222 Fruit trees and berry plantations Semi-open area 3 0
223 Olive groves Semi-open area 3 0
231 Pastures Agricultural area 2 0
241 Annual crops associated with permanent crops Agricultural area 2 0
242 Complex cultivation patterns Agricultural area 2 0
243 Land principally occupied by agriculture with significant areas of natural vegetation Agricultural area 2 0
244 Agro-forestry areas Forests 4 0
311 Broad-leaved forest Forests 4 0
312 Coniferous forest Forests 4 0
313 Mixed forest Forests 4 0
321 Natural grasslands Semi-open area 3 0
322 Moors and heathland Semi-open area 3 0
323 Sclerophyllous vegetation Semi-open area 3 0
324 Transitional woodland-shrub Semi-open area 3 0
331 Beaches dunes sands Semi-open area 3 0
332 Bare rocks Semi-open area 3 0
333 Sparsely vegetated areas Semi-open area 3 0
334 Burnt areas Semi-open area 3 0
335 Glaciers and perpetual snow Semi-open area 3 0
411 Inland marshes Wetlands 5 0
412 Peat bogs Wetlands 5 0
421 Salt marshes Wetlands 5 0
422 Salines Wetlands 5 0
423 Intertidal flats Wetlands 5 0
511 Water courses Water bodies 6 0
512 Water bodies Water bodies 6 0
521 Coastal lagoons Water bodies 6 0
522 Estuaries Water bodies 6 0
523 Sea and ocean Water bodies 6 0
999 NODATA NODATA 9 0

With a sort-of ‘raster attribute table’ in place, we can now reclassify the detailed landuse classes to the broader categories listed above.

This leaves us with the following count of ~100x100m cells per land use category

Herwijnen Den Helder
Urban area 1860678 809932
Agricultural area 6256242 3384907
Semi-open area 194271 133709
Forests 1437047 498510
Wetlands 346616 364301
Water bodies 3851726 8890353

3.2.4 Resampling the land use map to a lower resolution

The cellsize of the PPIs is 500, 500x500, 500 meters, but the land use map is much more finely detailed (~100x100m), so we need to resample the latter to derive a land use map at a 500, 500x500, 500 meter resolution as well.

As the resolution of the PPIs is so much higher than that of the land use map, we need to resample the latter to a lower resolution. Instead of classifying a single pixel in the PPI as belonging to a single land use class, we will do so using land use proportions. We therefore calculate the proportions belonging to each of the land use classes within an area of cells roughly the size of the PPI pixels. Subsequently we resample this to match 1:1 with the PPI pixels and store the land use proportions for each of the land use classes in the PPIs. This is done in the calculate_land_use_proportion() function below.

We can now calculate the proportions and will save these to some raster files for potential inspection in GIS software.

By now the resampled land use raster should be very similar to the PPI raster, with the exception of — of course — the values contained within.

## [1] TRUE
## [1] TRUE

Ok, let’s save a copy of what we have so far.

3.4 Calculate distance to nearest urban area

We can use the distance to the nearest inhabited urban area as a proxy for disturbance. To calculate this, we reclassify the raster to cells containing inhabited urban area and everything else. For every cell on the raster that is not a cell we have just classified as inhabited urban area, we will calculate the distance (in meters) to the nearest cell classified as inhabited urban area.

And we add these values to the PPIs again.

3.5 Add population density

Another proxy for disturbance is simply the number of humans living in a certain area. The Dutch Central Bureau of Statistics (CBS) annually publishes a dataset containing the number of inhabitants organized in 500x500m grid cells. We will now add this to the PPIs as well.

## Reading layer `CBS_VK500_2018_v1' from data source `/mnt/volume_ams3_01/raw/population-density/2019-CBS_VK500_2018_v1/CBS_VK500_2018_v1.shp' using driver `ESRI Shapefile'
## Simple feature collection with 151108 features and 31 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 13000 ymin: 306500 xmax: 278500 ymax: 619500
## proj4string:    +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs
## Warning in .local(x, ...): all fact(s) were 1, nothing to aggregate
## Warning in .local(x, ...): all fact(s) were 1, nothing to aggregate

Once again, verify if the PPI pixels match up exactly with the CBS population grids

## [1] TRUE
## [1] TRUE

Now as a final step we will calculate the total population within the neighborhood surrounding the PPI pixels, to get a more representative measure of disturbance potential in the surrounding area.

And add it to the PPIs.

And finally we save these PPIs again.

3.5.1 Plotting human parameters on an interactive map

We have now generated ‘human relevant’ parameters for later modelling stages. Let’s plot them on an interactive map again for reference.

3.5.1.1 Interactive map of Herwijnen

3.5.1.2 Interactive map of Den Helder